t=0:pi/20:2*pi;
x1=sin(t);y1=cos(t);
plot(x1,y1,'r')
axis equal
h=0.08;x=-2:h:2;y=-2:h:2;
[X,Y]=meshgrid(x,y);
if X^2+Y^2>1
    Z=1./sqrt(X.^2+Y.^2);
else
    Z=-1./sqrt(X.^2+Y.^2);
end
[PX,PY]=gradient(-Z,h);
hold on
contour(x,y,Z,[-1.9,-1.2,-0.8,-0,5,-0.1,0.1,0.5,0,8,1.24,1.9],'b');
quiver(X,Y,PX,PY,'k')
